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Abstract 

The algorithm which generates the ATI time scale at the National Institute of Stan¬ 
dards and Technology (NIST) has generated a scale with many desirable properties 
since 1968. Five of these are as follows: 

1. The fractional frequency variation of the scale is smaller than any clock in the 
scale for all integration times. 

2. The algorithm adaptively estimates the weights of clocks in real time. 

3. The scale is much more reliable than any individual clock. 

4. One can add or remove clocks from the scale easily, with a minimum impact on 
the scale. 

5. One can correct the ensemble for calibrations against a primary reference. 

There are three other properties we would like to obtain: 

1. Automatic frequency step detection. 

2. A scale optimized for post—processing, including running both forwards and back¬ 
wards in time. 

3. A scale that can run with minimal supervision for use in non-technical environ¬ 
ments. 

It turns out that simply estimating a variance of the frequency state of the clocks 
facilitates all three of these new properties. We report here a new algorithm which uses 
techniques from Kalman filtering to estimate this variance. Results from simulation 
and applications to real clocks arc presented also. 


INTRODUCTION 


Ideally, a time scale algorithm samples an ensemble of clocks to generate time arid frequency with 
more reliability, stability, and frequency accuracy than any of the individual clocks in the ensemble. 
In this paper we study an approach to this ideal. 

A time scale algorithm calculates the time offset of each of the ensemble clocks at a given reference 
time. Ensemble time, the time of the scale, is realized by applying the appropriate correction to any 
one clock. If there is no measurement noise this value is independent of which clock is used. The input 
to the algorithm at a given reference time is the time difference between each clock and a particular 
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clock. The algorithm also requires estimates of the deterministic arid stochastic parameters which 
characterize each clock’s frequency offsets. 

It is important to notice that the time of a clock is not measured. One measures only time differences 
between clocks. Thus the ensemble time which the algorithm generates is not observable. For this 
reason it is inappropriate to use an accuracy algorithm, such as a Kalman filter, to generate time 
by minimizing time error. We can, and do in the algorithms discussed here, optimize time and time 
interval stability. 

It is also important to realize that a clock as a physical system produces a frequency. The time of a 
clock is artificially derived from the frequency, which is the true physical quantity. Because of this, all 
the parameters which characterize clock performance describe aspects of the frequency. One can use 
these parameters to optimize time uniformity and frequency accuracy. An algorithm that optimizes 
time accuracy should simply allow the clock with the best long term stability to dominate the scale, 
thus sacrificing much of the performance of other clocks, especially in short term. We will look at 
these things further as we go along. 

The ATI time scale algorithm at the National Institute of Standards and Technology (NIST) has 
generated a scale with many desirable properties since 1967. Five of these are as follows: 

1. The fractional frequency variation of the scale generally appears smaller than any clock in the 
scale for all integration times. 

2. The algorithm adaptively estimates the weights of clocks in real time. 

3. The scale is much more reliable than any individual clock. 

4. We can add or remove clocks from the scale easily, with a minimum impact on the scale. 

5. We can correct the ensemble for calibrations against a primary reference. 

There are three other properties we would like to obtain: 

1. Automatic frequency step detection. 

2. A scale optimized for post-processing, including running both forwards and backwards in time. 

3. A scale that can run with minimal supervision for use in non-tcchnical environments. 

The new algorithm we report here combines aspects of the NIST ATI algorithm with techniques from 
Kalman filtering to estimate clock states representing the random walk plus drift of the frequency 
offsets of each clock, as well as the variance of these states. For a given clock, this state, which we call 
“Y”, is not a physical state, but a mathematical estimate of the frequency offset in the presence of 
white noise modulating the frequency. The variance of this state gives us a confidence of this estimate. 
Having this estimate facilitates attainment of all three of our goals. 

Frequency step detection always requires examining the data over some time period. In nearly real 
time operation it is possible to compare the estimate of average frequency offset over an interval with 
the filtered estimate from the beginning of that interval. Using the estimate of frequency deviation as 
a test for outliers, we can determine if a frequency step occurred in the recent past. 


344 



We can also smooth our estimates of “Y” in post-processing by combining the forward and backward 
filters. The proper way to combine these is well-defined in Kalman filter theory. Essentially, we use 
the reciprocal of the forward and backward filter variances as weights to combine their respective 
state estimates at a given reference time. We must be careful not to incorporate the data at that 
time in both filters before combining them, else they will not be independent estimates. Thus, at a 
given reference time, we use the extrapolated estimates of state arid variance from one direction, the 
backward filter for example, and combines this state estimate with the one from the forward direction 
which has been updated with the data. 

The ATI time scale is adaptive once we know the clock parameters characterizing the short and long 
term behavior of the clocks. Since our new scale uses the Kalman formalism, wc can enter a new clock 
and the algorithm will optimally adapt its estimate of frequency offset variance. This allows us to 
enter new clocks without perturbing the scale. 

Both results from simulation and applications to real clocks are presented. 


THEORY: ATI 

We first present the ATI algorithm and then show how we have modified it. In the ATI algorithm each 
clock has two states which are estimated: the time and frequency offsets of the clock from ensemble 
time. Frequency drift can be entered and used, but it is not estimated adaptively by the algorithm. 
The ATI algorithm is a three-tiered process: the time update, the update of the variance of the time 
offset, and the frequency update. The weight of a clock is proportional to the reciprocal of the variance. 
Each of these updates can be broken down into two steps: an initial estimate and the update. 


1. Time Prediction: 


Xi(t + r) = X,{t) + (Y4i)t Ar/2.r (l) 

We predict the time offset from ensemble time, X ,, of clock i for the current measurement time (t + r) 
based on the previous estimates at time t of time offset, A 7 ,-, and filtered frequency, Yj, and the entered 
frequency drift, Z?,-. 

2. Time Update: 


X^t + r) - £ = -t- r) - X {j (l - (2) 

We update the time offset of each clock j against the scale at time t f- t given the measurements 
Xij(t -1- r). Clock j 's offset is estimated using the measurement and prediction of each other clock, i, 
then these estimates are combined in a weighted average. The weights ?/;,■ are determined adaptively 
in equations 6 - 10. This is the maximum likelihood estimate of Xj if the are proportional 
to the reciprocal of the variances of the time residuals, and the residuals have a gaussian normal 
distribution.! 1 ! 
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3. Frequency Estimate: 


fi(t + r) = 


Xi(t + r) 


Xj(t) 


(3) 


Yi is the estimate at time t + t of the average frequency of clock z over the interval t based on the 
latest two time updates, X{. 


4. Frequency Update: 


U(< + 0 = + 0 + m i y i( t )l W 

We incorporate the previous frequency update into an exponentially filtered estimate of the current 
average frequency offset of clock t. The exponential frequency-weighting time constant (m,) is deter¬ 
mined from the relative levels of white noise FM and random walk (or flicker) FM for clock i (equation 

5)- 


5. Frequency Update, Exponential Time Constant: 


1 / 2-1 + 



(5) 


We determine ra, used in equation 4 to form the filtered estimate of the frequency of clock z. Here, 
T Mini is the integration time which gives the minimum value on a cr y {r) plot given that the clock’s 
stochastic deviations are characterized by white and random walk FM. to is the minimum t value 
used for computing This value of zrq can be shown to optimize the stability in predicting time 

(equation l) given these two kinds of noise in the clock (white and random walk FM)W. If white 
FM and flicker FM are more suitable models, then m,- can be approximated as t x /tq, where r x is the 
intercept value of r on a (7 y (r) plot for the white and flicker FM. 


6. Variance Estimate: 


€i [fl t + r) - Xi (t + r) I + Ki (6) 

£,• is the lack of predictability of clock i over the interval r, being the difference between the prediction 
and the update. Thus it is an estimate of the deviation (square root of variance) of the clock based on 
the current measurement cycle. The additive term Ki accounts for the fact that the term in brackets 
on the right-hand side of equation 6 is biased because clock z is part of the ensemble. See equation 
10 to calculate if,-. 
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7. Variance Update: 


e?(r) estimates the mean squared time error of clock i by filtering exponentially t 
*’s deviation from the current measurement cycle. Since the noise characteristics 
be stationary, past measurements are de-weighted in the filtering process. The ti 
filter is typically chosen to be N T = 20 days for cesium clocks, representing the t 
white FM level to be constant. The initial value of e?(r) can be estimated as r 2 a 


xponentially the estimate of clocl 


8. Ensemble Variance: 


£ * (r) Is IS) I 

cl(r) forms an estimate of ensemble time error, Any clock can 
performing clock cannot harm the stability of the ensemble. 


9. Adaptive Clock Weights; 


W{ is the weight to be used in equation 2 for clock When calculated this way, the resulting ensemble 
time stability can be shown to be optimized in a maximum likelihood sense, assuming a normal 
distribution of the noise of clock i with variance e?(r). 

10. Bias of the Error Estimate: 


Ki estimates the bias in the error estimate from the firs 
estimate is biased small, on the average, because each do 
through its weighting factor. The larger a clock’s weight, 
of a normal distribution of clock noise the size of the bi 
which is added to equation 5 in order to remove the bias 


term ori the right of equation 6. This error 
k. is a member of the ensemble and sees itself 
./he larger is the bias. Under the assumption 
5 can be estimated as given by equation 10. 


THEORY: AT1-PLUS-VARIANCE 


What is missing here is an estimate of the variance of the resic 
ensemble, Y. In our approach, we interpret .X, as a measuremi 










scale. Thus the first difference, Y|- of X{ 7 as in equation 3, is a measurement of the frequency offset of 
clock i from the scale. We use a simple Kalman formalism to filter this measurement to estimate Y, 
where we use e t (r) as noise of the measurement. We model Y as having a random-walk noise and a 
fixed drift. Thus, Y is not the physical frequency as produced by the clock, but since we are filtering 
down the white frequency modulation (FM), it is only the random-walk component of the frequency 
of the clock, plus any drift. 

The X terms do reflect the physical time offset of the clock from the scale, thus incorporating the 
white FM. The Y term is used to better predict the X values. We substitute equations 4 and 5 with 
the Kalman equations: 

11. System Model: 


Y,.(t + r) =Yi(t) + Ar + »?(r). (11) 

Y{ is the random walk component of the frequency offset of clock i from the scale plus the drift offset. 
The random walk is driven by the white noise process r}(r). 

12. Measurement Model: 


£■ = *;■ + «(»■) (12) 

The measurement Yi is a direct measurement of Y,, plus white noise. 

The actual equations used for update are: 

13. Variance Prediction: 


P(t + r) = P(i) + ** * r. (13) 

P is the prediction of variance of the residuals of Y which grows with r according to cr^, the variance 
of the white noise process, q, driving the random-walk FM. 

14. Frequency Update: 


T(t + r) = 


4(7-)*y(t) + P*K 

e?(r) + P 


(14) 


We see that the Kalman formalism also gives us an exponential filter on Y. Thus in steady state this 
algorithm reduces to ATI if the weights are chosen properly. This implies that this algorithm inherits 
the ability that ATI has to model flicker frequency. 
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15. Variance Update: 


P- 


<?(0 * P 

e t 2 (r) + P 


(15) 


In the Kalman formalism the system parameters are known in advance. This system is a modification, 
an adaptive Kalman filter, where we estimate the “measurement noise,” e‘-(r). This allows the variance 
of the residuals of Vto evolve both from an initial value, as is normal for the Kalman filter, and if 
ef(r) changes. This allows the exponential filter parameters on Y, as expressed in 14, to change with 
time, both after entering a new clock, and if the white FM level of the clock changes. 


It is possible to solve for the steady state form of these equations and make identifications between 
the ATI algorithm and our new algorithm. We find that the steady state value of P is: 


16. 


P = 


fM 

v 2 


(16) 


Making the appropriate identifications between equations 4 and 14 we find: 


17. 

m 


P -1- 7? 2 


(17) 


Using 16 in addition yields: 


, , e 2 (r) 

m(m + 1) = —— 

1 r 


(18) 


These equations allow us to compare the performance of the two algorithms with the parameters m 
and r] set consistently. 

We close this theoretical section mentioning results from elsewhere. Jones and Tryout have designed 
a time scale algorithm which is purely a Kalman filter. This scale, called TA(NIST), has been run 
at NIST in parallel with ATI since about 1983. That filter is mathematically identical with the 
ATI algorithm for the time and frequency predictions and updates^ 3 !. The difference among these 
algorithms is the weighting of clocks in the time update and the exponential filter parameters in the 
frequency update. These differences effect the ensemble time they generate, which is realized as the 
time offsets, X,-, of the clocks against ensemble time. 

We will show in simulation and with real clocks that the pure Kalman filter time scale sacrifices 
short term performance, and simply follows the clock with the best long term performance. This is 
consistent with the design of Kalman filters in general which minimize error. The Jones Tryon filter 
attempts to estimate and minimize both time and frequency error. An additional problem with this is 
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that since time is unobservable, elements of the covariance matrix grow without bound. In practice, 
with a good ensemble of clocks, this growth is not large enough to cause computer overflow errors in 
any reasonable amount of time, though it is suggestive of an undesirable situation. 

SIMULATION 

In simulation we show the following: 

1. Both algorithms ATI and ATI plus frequency variance produce a time scale apparently better 
than the best clock in the scale at all integration times. 

2. The TA(NIST) algorithm is dominated by the clock with the best long term performance at all 
integration times. 

3. The ATl-plus-frequency™ variance estimate of the confidence on the frequency offset estimate 
appears to be a reasonable estimate. 

4. The use of this confidence estimate to determine frequency steps improves long term performance 
of the time scale. 

Figure 1 illustrates item (l). Here we have generated data simulating clocks with various levels of 
white FM and random- walk FM. We have treated the problem as we would with real clocks where we 
only measure clock differences. We have computed the stability of each clock using an N-cornered hat 
technique^ 4 !. The stability of the scale we have determined by taking the output value of clock minus 
scale and subtracting the generated value of clock minus truth. If we look directly at the variances of 
the generated data, we see significant differences between the variances computed directly, and those 
estimated from N-cornered hat (Figure 2). These differences must be due to apparent correlations 
in the data. This should come either from the finite data length, or from real correlations in the 
pseudo-random number generator. If the generated clocks are truly correlated, then the algorithm 
can only produce a variance better than the uncorrelated part. Wc notice that the scale seems to 
follow the shape of the variances from the N-corner hat. This suggests correlation in the generated 
data. 

Figure 3 shows a comparison of the output of a version of the Kalman filter which defines TA(NIST) 
with the simulated input data. We see that the Kalman algorithm has t,hc stability of the clock with 
the best long term variations. 

Figures 4 and 5 show the residuals from the ATl-plus-variance algorithm compared with the estimated 
confidence, from the estimated variance of frequency residuals. The algorithm estimates the random- 
walk component of a clock’s frequency offset from the ensemble time. Since these are generated clocks 
we know the true value of the random-walk component of frequency of that clock versus the true value 
of the time scale. The differences of these two, the estimate minus truth, arc the residuals plotted. 
The sigma value used in the plot is the root-mean-square of the estimated deviation of the clock plus 
the estimated deviation of the scale. The line plotted is the three-sigma value. This should be a 99.8 
percentile. Over the 700 points plotted we should get 1 or 2 residuals crossing the lines of the sigmas. 
This seems to be the case. 

Last, we inserted frequency steps in the simulated clocks. Figure 6 shows the frequency offset from 
the scale of the simulated clock 9, with a frequency step of 1 x 10 -12 on MJD 46500. This clock 
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was given a white FM level of 30 ns, and a random walk FM level of 0.5 ns, both at 1 d. Figure 7 
shows the estimate from the ATl-plus-variance scale of the random-walk component of frequency. 
The reduction in the white FM is apparent. The scale was able to detect the frequency step, with the 
step detector set at 4 sigma, with sigma defined as above. When such a step is detected, we re-run 
the scale, removing the clock with the step until the scale can learn the new frequency value. Figure 8 
shows the frequency offset from the scale of the simulated clock 1, with a frequency step of 2 X 10 “ 12 
on MJD 46100. The noise of this clock is almost all random walk FM as compared to its level of white 
FM. The estimate from the scale shows very little smoothing. Yet, even in this case, the frequency 
step detector automatically found the step and removed the clock from the scale. In Figure 9 we 
see the benefit from having detected the frequency steps. There is a significant improvement at an 
integration time of 128 d and longer. 


REAL DATA 

Last, we took data from real clocks at NIST over the period from December 31, 1988, to October 
30, 1989. We ran our ATl-plus-variance algorithm on this data, including automatic frequency step 
detection and recalculation. We also took data using GPS common view measurements! 5 ! with other 
laboratories: PTB, USNO, TUG, and NRC. Using an N-cornered hat technique, we were able to 
determine the variance of each of these. We compared these results with a similar analysis using the 
official NIST ATI time scale. The results for the two scales run on NIST clocks as well as USNO 
and PTB are plotted in Figure 10. We find that the official ATI and the new ATl-plus-variance 
scale are similar, though in long term the official scale is somewhat better. The official ATI scale is 
watched carefully and administratively checked for time and frequency steps, as well as changes in 
clock performance in general. We find that human care adds much to a time scale. 
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ATI + Variance on Simulated Clocks 



Figure 1: We have computed the stability of each simulated clock using an N- 
cornered hat technique. The stability of the scale we have determined by 
taking the output value of clock minus scale and subtracting the generated 
value of clock minus truth. The ATI scale outperforms all clocks at all 
integration times. 

Truth vs N-Comer Hat 





• 9 

S3 


’*• 10* Ilf 10' 

Figure 2: If we look directly at the variances of the generated data, we see 
significant differences between the variances computed directly, and those 
estimated from N-cornered hat. These differences must be due to apparent 
correlations in the data. 


Kalman on Simulated Clocks 



10 7 10* 


Figure 3: A comparison of the output of a version of the Kalman filter which 
defines TA(NIST) with the simulated input data. We see that the Kalman 
algorithm has the stability of the clock with the best long term variations. 
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Simulation Clk 1: Residuals vs. Estimated Confidence 

White FM = 1 ns. Random Walk FM = 15 ns. both @ 1 day 



46000 46100 46200 46300 46400 46500 46600 46700 

Time, MID 


Figure 4: The residuals from the ATI-plus-variance algorithm compared with 
the estimated confidence, from the estimated variance of frequency residuals. 
The line plotted is the three-sigma value. 


Simulation Clk 9: Residuals vs. Estimated Confidence 

White FM = 30 ns, Random Walk FM = 0.5 ns, both @ 1 day 



46000 46100 46200 46300 46400 46500 46600 46700 


Figure 5: The residuals from the ATI-plus-variance algorithm compared with 
the estimated confidence, from the estimated variance of frequency residuals 
The line plotted is the three-sigma value. 
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Frequency-Step Detection Improves Long-Term Stability 
Simulated Clocks 



Time, seconds 


Figure 6: The frequency offset from the scale of the simulated clock 
a frequency step of 1*10“ 12 on MJD 46500. 


Comparison of Ensembles 



Time, seconds 


re 7: Th 


late from the ATI-plus-variance scale of the rando 






Figure 8: The frequency offset from the scale of the simulated clock 1, with a 
frequency step of 2*10" 12 on MJD 46100. The estimate from the scale shows 
very little smoothing, yet the frequency step detector automatically found the 
step and removed the clock from the scale. 



Figure 9: The benefit from having detected the frequency steps is seen in 
this diagram. There is a significant improvement at an integration time of 
128 d and longer. 



Figure 10: The results for the two scales, the official ATI and the new ATI- 
plus -variance , run on NIST clocks as well as USNO and PTB. 
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QUESTIONS AND ANSWERS 

GERNOT WINKLER, USNO: There are quite a few things that I didn’t understand, but there 
is one thing that bothers me—at the beginning you said that, since you don’t measure your time scale, 
but generate it, the use of a Kalman filter is inappropriate. Then you turn around and use it. Could you 
elaborate on that? 

MR. WEISS: The use of a pure Kalman filter is inappropriate. Using Kalman techniques in part of 
the time scale is very useful. This is not a criticism of Sam Stein’s algorithm which does a very similar 
thing. Kalman techniques have very important places. But to use the Kalman to generate the time itself 
is inappropriate. The result of that is what I have shown. You lock onto the clock with the best long term 
stability and you lose the short term. 

MR. STEIN: I think that I would like to make a comment on that myself. The Kalman filter is merely 
a computational technique. It doesn’t define a time scale, in the sense that we use the word. Whether or 
not you can use a Kalman filter, I strongly disagree with your statement, is simply whether or not you can 
formulate an appropriate definition of the time scale, definition of the ensemble in a way that is compatible 
with requirements. In fact, troubles that you pointed out, which are real, in some of the past Kalman filters 
were a result of the incorrect identification of the system noise covariance. An identification that did not 
agree with what the people that did those time scales really wanted to get out of the system. Other noise 
models produce different results. 

JAMES BARNES, How often do you see frequency steps? 

MR. WEISS: Dave Allan can probably answers that. 

MR. ALLAN: It is a fairly difficult question because we don’t really understand the mechanism. We 
do know that they occur. I have seen some clocks where they occur maybe once a year, several times a year 
and some clocks that seem not to have step problems. They can be fairly large and fairly significant. 

MR. WEISS: We have recently discovered a possible source of the steps which was humidity, and we 
had never realized that before. Tracking and controlling humidity may remove some steps that we didn’t 
understand the cause of before. As we do more research, we may understand the causes of other frequency 
steps. 

UNIDENTIFIED QUESTIONER: How do you define ‘truth’? 

MR. WEISS: I don’t define truth. When I generate a clock, I know what I have put into it. I generate 
a clock as being an offset against zero, zero being truth. If I make it up, I know what truth is. 

SAME PERSON: Is it related to the Planck time? 

MR. STEIN: No. You can relate it to what Einstein said, ‘Time is what a clock reads’. 

UNIDENTIFIED QUESTIONER: I have a question for you, Sam. You said that your model 
differed form the others in the process noise. As far as I can tell, the problem is one of observability. You 
can’t change the observability of the model by changing process noise. If your H matrix is the same as 
everyone else’s, then you have the same observability problem they do and all you are doing is de-weighting 
somehow past data. I think that I know how you are doing that and I can tell you. 

MR. STEIN: The problem, as Mark pointed out, is not a problem of observability. We can discuss it 
further off-line. 

DONALD PERCIVAL, U OF WASHINGTON: I just had a question to clarify what you are 
doing here. You are not using the (indecipherable)-Jones Kalman filtering model, is that correct? What 
you are doing then is taking one component of that model, the random walk component, and then using a 
Kalman filter by itself, a very simple one, to estimate that. 

MR. WEISS: Exactly. I am using the output of ATI time estimate against the scale as if it were a 
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measurement of frequency against the scale. Then I filter that to obtain only the random walk component 
plus the variance, or confidence, on that estimate. 

MR. PERCIVAL: Is there any problem, do you feel, with the observability. For example, in the time 
model, the random walk component is only one component of what the thing is actually doing. To separate 
that out, you are saying that it can be filtered out in a certain fashion. 

MR. WEISS: I think that that is a subtle question, because the. scale .is a mathematical construct. 
It is something that, mathematically you are making up. The only observability of that is within the scale, 
there is no physical measurement of the clock against the scale. The closest; that we can get to it is what 
the scale says is the estimate. There is no question of observability because T am looking at something that 
is defined. 

GERNOT WINKLER, USNO: I would like to go back to the frequency step detection thing. It 
seems to me that it is intrinsically related to the question of post-processing or not;. You can readily make 
a decision as to whether or not there was a frequency step, after the fact. You define as a frequency step 
as a systematic change in the average frequency. Whether it is systematic or random can only be decided 
with any kind of confidence only after the fact for small frequency steps. It is a question of the size of the 
frequency steps that you want to estimate vs. the question of post-processing or immediate detection. Here 
you have to solve the problem of whether you can tolerate post-processing or not. That will determine 
whether you can reliably detect frequency steps and account for them. 

MR. WEISS: That’s true. But, if you can detect that a step happened, real time, that is, you 
detect now that a step happened yesterdaj—and the step is large enough, can get some benefit from 
that now, in real time. You get the most benefit in post-processing, because then you go back over the last 
day as well. If the step is large enough that it is pulling the scale now, yon can readjust that. You can 
also re-compute the scale for the last day and, even if you can’t use that in real time, you can say that the 
frequency of the scale is now off, somewhat, from what I am actually outputting, and steer back to what the 
frequency should be. 
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